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ABSTRACT 



A finite element formulation was used to solve the steady- 
state two-dimensional conduction heat transfer equation in 
the condenser wall section of an internally finned rotating 
heat pipe. A FORTRAN program using this method was coupled 
with the COPES/CONMIN program for sensitivity analysis of 
design variables and for automated design of the internal 
heat pipe geometry. 

With water as the working fluid, numerical results ob- 
tained for copper and stainless steel heat pipe condenser 
sections indicated that for the maximum heat transfer rate, 
the designer should machine as many fins as the condenser 
material and the manufacturing process will allow. A saw 
tooth profile is preferable to spacing between fins. 
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I. 



INTRODUCTION 



A. THE ROTATING HEAT PIPE 

The rotating heat pipe is a closed container designed 
to transfer a large amount of heat in rotating machinery. 

Its three main component parts are: a cylindrical evapora- 

tor, a truncated cone condenser, and a working fluid as 
shown in Figure 1. 

At rotation above the critical speed of a rotating heat 
pipe, the working fluid forms an annulus in the evaporator, 
and will be vaporized by heat addition to it. The vapor 
flows toward the condenser as a result of a pressure differ- 
ence, transporting the latent heat of vaporization with it. 
External cooling of the condenser causes the vapor to con- 
dense on the inner wall and release its latent heat of eva- 
poration. The centrifugal force due to the rotation has a 
component acting along the condenser wall that will act to 
drive the condensate back to the evaporator where the cycle 
is repeated. 

In a conventional heat pipe, the force driving the con- 
censate back to the evaporator is due to capillary action, 
which poses a limit to its operation. The rotating heat 
pipe is not limited by capillary action and, unlike the 
thermosyphon which depends on gravity to cause condensate 
return, can be used in any orientation [1]. 
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B. OPERATING LIMITS OF A ROTATING HEAT PIPE 



The first theoretical investigation of the rotating 
heat pipe conducted at the Naval Postgraduate School was per- 
formed by Ballback [2] in 1969. He studied the limitations 
of performance imposed on the rotating heat pipe due to 
various fluid dynamic mechanisms. Using existing theory and 
experimental correlations, he was able to estimate the sonic 
limit, boiling limit, entrainment limit, and the condensing 
limit of performance. 

1. The Sonic Limit 



pipe, it is possible to reach a limiting flow rate of the 
vapor brought on by a choked flow condition in the pipe. 

This condition imposes a limiting value on the amount of 
energy the vapor can transport, thus reducing the effective- 
ness of the heat pipe. The limiting heat transfer rate 
becomes 



When increasing the heat flux in a rotating heat 




( 1 ) 



and the vapor velocity is considered to be sonic, 



U v = c = / g 0 kRT 



( 2 ) 



where 



U v = velocity of the vapor in ft/sec, and 

2 

A = cross sectional area for the vapor flow in ft 
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c = sonic velocity in ft/sec 

2 

gg = gravitational constant, 32.1739 f t-lbm/lbf-sec 
k = ratio of specific heats 
R = gas constant in ft-lbf/lbm ?.R, and 
T = absolute temperature in ?.R . 

2 . Boiling Limit 

Kutateladze [3] postulated that the transition from 
nucleate to film boiling is totally a hydrodynamic process. 
He determined a theoretical formula for predicting the burn- 
out flux 

Q t = K /p^ A b h fg {a g(p f - P v )} ±/4 (3) 

where 

K = constant value 

. 3 

P v = density of the vapor in lbm/ft 

2 

A b = heat transfer area in the boiler in ft 
h^ = latent heat of vaporization in Btu/lbm 
a = surface tension in lbf/ft 

2 

g = acceleration of gravity m ft/hr 

3 

p p = density of fluid in lbm/ft 

3 

P v = density of vapor in lbm/ft 
The experimental data obtained by Kutateladze suggested a 
value for K in the range of 0-. 13 to 0.19. 
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3 . Entrainment Limit 



The flooding constraint in a wickless heat pipe was 
examined by Sakhuja [4] who developed the correlation 



<4 



A x g2 h fz / g^Pf-Py’Py 

c . , * / -v 1 / 4 t 2 

{1 + (p v /p f ) } 



(4) 



where 

= heat transfer rate in Btu/hr 

2 

A = flow area m ft 
x 

C = dimensionless constant, 0.725 for tube with sharp 
edged flange 

hjr = latent heat of vaporization in Btu/lbm 

2 

g = acceleration due to gravity in ft/hr 

D = inside diameter of heat pipe in ft 

3 

= density of the fluid in lbm/ft 

3 

P v = density of the vapor in lbm/ft 

4 . Condensing Limit 

Ballback [2] determined the condensation solution 
for a rotating heat pipe by modeling the condenser section 
of a rotating heat pipe as a rotating truncated cone. He 
developed the following expression for the condensation 
limit : 




P f u 2 h fK {T s -T w } 3 

. 2 
sm <p 



1 

4 



} 



1 11 

{[R + Lsin d> ] 3 -R 3 } 4 

o Y o 

( 5 ) 
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where 



Q 

k 

P 



t 

f 

f 



w 




<f> 

R 

L 



y 



f 



= total heat transfer rate in Btu/hr 

= thermal conductivity of the condensate film 
in Btu/hr-ft-°F 

3 

= density of fluid in lbm/f t 
= angular velocity in 1/hr 
= latent heat of vaporization in 3tu/lbm 
= saturation temperature in °F 
= inside wall temperature in ?F 
= viscosity of fluid in lbm/ft-hr 
= half cone angle in degrees 
= minimum wall radius in ft 

= length along the wall of the condenser in ft 
= viscosity of the fluid in lbm/ft-sec 



The condensing limit equation is a function of the geometry 
and speed of the rotating heat pipe, and the physical pro- 
perties of the working fluid. 

Tantrakul [5] calculated these limitations for a heat 
pipe with specific physical characteristics as shown in 
Table I, with the results shown in Figure 2. 



TABLE I 



Specification of a Typical 


Rotating 


Heat Pipe 


Length 


14.000 


inches 


Minimum diameter 


2.000 


inches 


Wall thickness 


0.125 


inches 


Internal half angle 


1.000 


degree 


Rotating speed 


2700 


RPM 
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Obviously from the results in Figure 2, the condensing limit 

t 

is the predominant limitation for the amount of heat that can 
be transferred from the heat pipe. However, the other 
limitations may become important as the heat pipe geometry and 
operating conditions are varied. 

In order to augment the heat transfer capacity of the 
heat pipe, recent efforts have been aimed at raising the con- 
densing limit line which may be accomplished by: 

a. a high value of cone angle, to increase the centrifugal 
driving force, 

b. some type of promoter of dropwise condensation to in- 
crease the value of the inside heat transfer coeffi- 
cient, h, or 

c. use of an internally finned condenser to increase 
the inner wall surface area and the value of h, since 
the presence of a fin will decrease the effective 
condensate film thickness. 

A high value of cone angle means a departure from the 
cylinder or shaft shape. Since the principal known appli- 
cation for the rotating heat pipe is in the cooling of rota- 
ting machinery this approach was not pursued. Although 
effective promoters of dropwise condensation exist, none as 
yet can be considered to be permanent, and this approach was 
likewise ruled out. The remaining alternative, using inter- 
nal fins in the condenser section to raise the condensing 
limit, was seen as the best choice. 
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C. ANALYSIS OF THE INTERNALLY FINNED ROTATING HEAT PIPE 

Pursuing the addition of internal fins as a way to raise 
the condensing limit, Schafer [6] developed an analytical 
’ model for a heat pipe with a triangular fin profile as shown 
in Figure 3. He assumed one-dimensional heat conduction 
through the wall and fin. Corley [7] for this same case 
developed a two-dimensional heat conduction model using a 
Finite Element Method, and also assumed a parabolic tempera- 
ture distribution along the fin surface. His results indi- 
cated a significant improvement in heat transfer performance 
of about 75% above that predicted by the one-dimensional model 
of Schafer [6]. However, Corley [7] cautiously noted a 
probable error of 50% existed at the fin apex, and conse- 
quently mentioned that there may be a total heat transfer 
error of as high as 15%. Tantrakul [5] modified Corley's 
computer program by increasing the number of finite elements 
in order to minimize the heat transfer error at the apex of 
the fin. His results with this modification converged with 
the results of Corley. Purnomo [1] developed a two-dimen- 
sional Finite Element Method solution using a linear trian- 
gular finite element model as shown in Figure 4. Purnomo ' s 
Cl] Finite Element Method program also worked and converged. 
Purnomo 's [1] code, when made to approach the geometry of a 
smooth tube, did not agree with the analytical and experi- 
mental data obtained by Schafer [6] for a smooth tube. This 
cast doubts about the validity of Purnomo ' s code. 
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The parametric studies conducted using Purnomo ' s code 
gave no clear indication of the best condenser geometry to 
maximize heat transfer. Also, his code was tedious to use 
and required numerous runs to obtain data since it was 
written to perform only one analysis at a time. 

Clearly a computer program that could make numerous runs 
with minimal data input and could also automatically find 
improved designs would be valuable. 



D. THESIS OBJECTIVES 

The objectives of this thesis were therefore: 

1. To modify Purnomo ' s [1] computer program so that it 
is compatible with the COPES/COMNIN program [8] and 
can be used for analysis and automated design of 
rotating heat pipes (internally finned or smooth). 

2. To compare results using Purnomo ' s code with analyti- 
cal results for a smooth tube obtained by Schafer [6] 
to determine if and where an error exists. 

3. To use the sensitivity analysis capability of COPES/ 
CONMIN in conjunction with the modified program to 
study heat transfer in an internally finned rotating 
heat pipe . 

4. To use the resulting program to obtain an optimum 
design for an internally finned rotating heat pipe 
to obtain experimental data to compare with the 
analytical results. 

5. To use the resulting program to obtain numerical 
results in place of data obtained from expensive 
experimental operations. 
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II. NUMERICAL OPTIMIZATION 



A. BACKGROUND 

Most design processes require the minimization or maxi- 
mization of some parameter which may be called the design 
objective. For the design to be acceptable, it must satisfy 
a variety of physical, aesthetic, economic and, on occasion, 
political limitations which are referred to here as design 
constraints. While part of the design problem may not be 
easily quantified, most of the design criteria can be 
described in numerical terms . 

To the extent that the problem can be stated in numerical 
terms, a computer program can be written to perform the 
necessary calculations. For this reason, computer analysis 
is commonplace in most engineering organizations. For 
example, in structural design the configuration, materials, 
and loads may be defined and a finite element analysis com- 
puter code is used to calculate stresses, deflections, and 
other response quantities of interest. If any of these 
parameters are not within prescribed bounds, the engineer 
may change the structural member sizes and rerun the program. 
The computer code therefore provides only the analysis of a 
proposed design, with the engineer making the actual design 
decisions. This approach to design, which may be called 
computer-aided design, is commonly used today. 
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Another common use of analysis codes is in tradeoff 
studies. For example, an aircraft trajectory analysis code 
may be run repetitively for several payloads , calculating 
the aircraft range, to determine the range-payload sensitivity. 

A logical extension to computer-aided design is fully 
automated design, where the computer makes the actual design 
decisions, or performs trade-off studies with a minimum of 
man-machine interaction [93. 

B. CONSTRAINED FUNCTION MINIMIZATION (CONMIN) 

Vanderplaats [10] developed an optimization program CONMIN 
capable of optimizing a very wide class of engineering pro- 
blems. CONMIN is a FORTRAN program, in subroutine form, that 
optimizes a multi-variable function subject to a set of in- 
equality constraints based on Zoutendijk's [11] method of 
feasible directions [12]. 

Three basic definitions are required to discuss the 
use of CONMIN: 

Design Variables - Those parameters which the optimization 
program is permitted to change in order to improve the 
design. Design variables appear only on the right hand 
side on an equation and are continuous. 

Design Constraints - Any parameter which must not exceed 
specified bounds for the design to be acceptable. 

Design constraints may be linear or nonlinear, implicit 
or explicit, but they must be continuous functions of 
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the design variables. Design constraints appear only 
on the left side of the equations. 

Objective Function - The parameter which is going to be 
minimized or maximized during the optimization process. 
The objective function may be linear or nonlinear, 
implicit or explicit, and must be a continuous function 
of the design variables. The objective function usually 
appears on the left side of an equation, but it may 
appear on the right side if it is also a design variable. 
Design constraints and objective functions are usually inter- 
changeable . 

C. CONTROL PROGRAM FOR ENGINEERING SYNTHESIS (COPES) 

Recall that the optimization program, CONMIN, was written 

in subroutine form. Vanderplaats [8] has developed a main 
program to simplify the use of CONMIN and to further aid in 
the optimization process. The user must supply an analysis 
subroutine named ANALIZ. What follows are programming guide- 
lines to ensure compatability with COPES. 

D. PROGRAMMING GUIDELINES 

In developing any computer code for engineering analysis, 
it is prudent to write the code in such a way that it is 
easily coupled to a general synthesis program such as COPES. 
Therefore, a general programming practice is outlined here 
which in no way inhibits the use of the computer program in 
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its traditional role as an analytic tool, but allows for 
simple adaption to COPES. This approach is considered good 
programming practice and provides considerable flexibility 
of design options. Only five basic rules must be followed: 

I. Write the code in subroutine form with the primary 
routine called as SUBROUTINE ANALIZ ( ICALC ) . The 
name ANALIZ is compatable with the COPES program and 
ICALC is a calculation control. Note that subroutine 
ANALIZ may call numerous other subroutines as required 
to perform the necessary calculations. 

II. Segment the program into INPUT, EXECUTION, and OUTPUT. 
The calculation control, ICALC, will determine the 
portion of the analysis code to be executed. ICALC=1; 
the program reads all data required to perform the 
analysis. Also, any initialization of constants 
which will be used repetively during execution is 
done here. This initial input information is printed 
here for later reference and for program debugging. 
ICALC=2; the program performs the execution phase of 
the analysis task. No data reading or printing is 
done here, except on user-defined scratch disc. Data 
may be printed here during program debugging, in 
which case it should be controlled by a print con- 
trol parameter which is read during input. In this 
way, this print may be turned off after the program 
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is debugged, but may be used again during future 
program expansion debugging. The reason that print- 
ing is not allowed during execution is that when 
optimization is being done, the code will be called 
many times with ICALC=2, resulting in voluminous 
print. ICALC=3; the results of the anlaysis are 
printed. Also the essential input parameters which 
may have been changed during optimization should be 
printed here for easy reference. In summary, when: 
ICALC = 1 Read input data. 

ICALC = 2 Execute the analysis. 

ICALC = 3 Print the results. 

III. Store all parameters which may be design variables, 
objective functions or constraints in a single 
labeled common block called GLOBCM. The order in 
which they are stored is arbitrary. A listing of the 
COPES program should be checked to see how many para- 
meters may be stored in GLOBCM (the dimension of 
ARRAY) . Initial distribution of COPES allows for 1500 
parameters . 

IV. During execution or output, no parameters which are 
read during input should be updated. For example, 
if variable X is initialized during input, the 
execution segment must not update X such as X=X + 

3.2. Instead a new variable, Y=X + 3.2 must be 
def ined . 
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V. Write all programs in standard language, avoiding 
machine dependent capabilities such as seven letter 
FORTRAN names (CDC). While this guideline is not 
essential to the use of the analysis code within the 
COPES program, it makes the analysis code much more 
transportable between different computer systems, a 
capability which easily justifies a slight reduction 
in efficiency on a given machine. 

Adherence to these guidelines not only leads to a more 
readable and machine independent computer code, but allows 
this code to be coupled to the COPES program without modifi- 
cation. 

Having written the analysis code, it may be executed 
either with a simple main program or within the COPES pro- 
gram to perform the analysis. To insure that guideline IV 
is followed, the following main test program is recommended. 
Note that this program calls ANALIZ twice with ICALC=2 and 
ICALC=3, to show that the same result is obtained repeti- 
tively . 

C MAIN PROGRAM TO CHECK SUBROUTINE ANALIZ. 

C READ, EXECUTE , AND PRINT 

DO 10 ICALC =1 , 3 
10 CALL ANALIZ ( ICALC ) 

C EXECUTE AND PRINT AGAIN TO BE SURE THE RESULTS 
C DO NOT CHANGE 
DO 20 ICALC= 2 , 3 
20 CALL ANALIZ (ICALC) 

STOP 

END 
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III. FINITE ELEMENT SOLUTION 



A. REVIEW OF THE PREVIOUS ANALYSIS 

Schafer [6] studied the one-dimensional model heat trans- 
fer solution and Corley [7] studied the two-dimensional 
model for an internally finned rotating heat pipe. Both 
used the same assumptions and boundary conditions based upon 
the analysis of 3allback [2], which are similar to those used 
in the Nusselt analysis of film condensation on a flat wall. 

The more important of those assumptions are: 

1. steady state operation} 

2. film condensation, as opposed to dropwise condensation, 

3. laminar flow of the condensate film along both the 
fin and the trough, 

4. static balance of forces within the condensate, 

5. one -dimensional conduction heat transfer through 
the film thickness (no convective heat transfer in 
the condensate film), 

6. no liquid - vapor interfacial shear forces, 

7. no condensate subcooling, 

8. zero heat flux boundary conditions on both sides of 
the wall section (symmetry conditions), as shown 

in Figure 5, 

9 . saturation temperature at the fin apex, 
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10. zero film thickness at the fin apex, and 

11. negligible curvature of the condenser wall. 

Purnomo [1] developed a two-dimensional Finite Element 

Method solution using a linear triangular finite element 
model as shown in Figure 4. Purnomo modif ied Corley's assump- 
tion that the fin apex was at the saturation temperature and 
allowed the value of the temperature at the apex to float. 

He assumed a parabolic temperature distribution along the 
fin surface. 

Purnomo ' s statement of the problem for the formulation 
of the Finite Element Method as shown in Figure 5 is 



2 2 

3 T x 3 T _ n 

2 2 “ U 
3x 3 y 



( 6 ) 



with the boundary conditions : 



a) along boundary 1 



-k ? -- h (T - T ) 
an 1 sat 



b) along boundary 3, -k P s h„ (T - T ) 

<3 n 2 00 



9 T 

c) along boundaries 2 and 4, r — = 0 

o n 



A detailed description of the numerical formulation is pre- 
sented in his thesis. 
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Purnomo ' s computer program consisted of a main program 
and three subroutines; 

a) the routine "COORD" used to define positions of 
system coordinate points, 

b) the routine "FORMAF" used to formulate the 
Finite Element Method equations, and 

c) the routine "BANDEC" as an equation solver for a 
symmetric matrix which has been transformed into 
banded form. 

B. THE COMPUTER PROGRAM DEVELOPMENT 

Purnomo ' s [1] two-dimensional finite element program 
is the basis for the present analysis code. The first 
task undertaken in the development of this thesis was to 
check Purnomo ' s code for validity. Theoretical results 
using Purnomo ' s code for the case of zero fins were compared 
to theoretical results and experimental data for a smooth 
tube obtained by Schafer [6]. Purnomo ' s results exceeded 
Schafer's results by a factor of approximately two. In 
studying this discrepancy, an error was discovered in 
Purnomo ' s code. In encoding formula (11-11) [1] for mass 
flow rate in FORTRAN, the' sin <f> ' term was dropped. The 
correct form of the equation is shown below. 

. p^u^(R n + xsin <|>)S*(x)2 sin <j> _ 

M tot (x) = [5 * (x) e+6 * ( X) tan a] 

3 y ^ 

( 7 ) 
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The effect of this error was that the film thickness along 
the length of the condenser remained very small. Purnomo ' s 
original code consequently indicated a much higher value of 
heat transfer rate than the correct analysis did. In a 
correct analysis code, the condensate film thickness grows 
continuously until the condensate reaches the evaporator. A 
plot of the condensate film thickness calculations from 
Purnomo ' s original code for the case of a smooth tube, his 
corrected code, and Schafer's program for a smooth tube are 
shown in Figure 6. Purnomo ' s corrected code for the case of 
zero fins agreed to within 8 percent of the results for 
heat transfer rate using Schafer's [6] analysis for a smooth 
tube . 

Since no experimental data existed for further comparison 
of the finned model, the next task undertaken was to adapt 
the analysis code to permit automated design and sensitivity 
analysis using COPES/CONMIN . Many modifications were made, 
some of which are mentioned here. The program was rewritten 
in subroutine form and segmented into input, execution, and 
output sections to make it compatible with COPES. Since 
COPES was written to use single precision mathematics and 
the analysis code uses double precision to allow for possible 
ill conditioning, the subroutine ANALIZ also makes the trans- 
formation from single to double precision. The initial value 
of the film thickness is calculated within the code based on 
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a formulation by Sparrow and Gregg [13]. Previously a 
constant value based on the same analysis was used. 

The modified code can be used alone for analysis of a 
given geometry or it can be used in conjunction with COPES/ 
CONMIN for a single analysis, sensitivity analysis, or 
automated design (optimization) . 

A listing of the revised computer program is included 
as Appendix B. 

C. SENSITIVITY ANALYSIS 

The sensitivity analysis feature of COPES/CONMIN was used 
to obtain data for the heat transfer rate Q as a function 
of various design variables. All design variables were held 
constant except the one of interest which was varied as 
specified in block Q of the COPES input data. 

For example, to obtain the heat transfer rate Q for 26 
different values of rotational speed, the following data 
input is required : 

$ Block P 

1 

15 

$ Block Q 

16 , 26 

3000., 100., 200., 300., 400., 500., 600., 700., 

800., 900., 1000., 2000., 3000., 4000., 5000., 6000., 

7000., 8000., 9000., 10000., 11000., 12000., 13000., 

14000., 15000., 3000. 
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For each plot only the variable on the ordinate is 
changed. All other design variables remain constant. The 
basic design used for sensitivity plots has the design 
variable values shown below: 



TABLE II 

BASIC DESIGN FOR SENSITIVITY ANALYSIS 



condenser length 


CLI 


= 


9 . 0 inches 


cone half angle 


CANGL 


Z 


1.0 DEGREES 


condenser radius 


RBASEI 


= 


0.775 inches 


wall thickness 


THICKI 


= 


0.03125 inches 


fin height 


B 


= 


0.025 inches 


speed of rotation 


RPM 


= 


3000 RPM 


saturation temperature 


TSAT 


Z 


100° F 


ambient temperature 


TINF 


= 


60° F 


outside heat transfer 
coefficient 


HINF 


z 


5000 Btu/HR FT 


fin half angle 


FANGL 


= 


10 


ratio of trough width 
to fin base width 


ZOA 


= 


12.8 


number of fins 


ZFIN 


z 


40.0 



Note: In the code, ZFIN appears on the left hand side of 

an equation and is therefore by strict definition not a 
"design variable" . ZFIN is calculated from the values of 
CBASE , EZERO, and EPSO. 

D. DESIGN OPTIMIZATION 

Counting parameters such as external heat transfer 
coefficient HINF, there are thirteen possible design 
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variables . Nine of these are geometric or functional para- 
meters such as wall thickness, fin height, and speed of 
rotation. The design variables, possible constraint func- 
tions, and the objective function appear in the Global 
common block GL0B1 in the code and are listed below by fortran 
name for clarity. 



DESIGN VARIABLES 
CLI 
CANGL 
RBASEI 
R2I 

THICKI 

3FIN 

TZ 

TSS 

TINE 

HINF 

FANGL 

ZOA 

RPM 

CONSTRAINT FUNCTIONS 
ZFIN 
BOA 
DIFF 

OBJECTIVE FUNCTION 
QTOT 



There are a wide variety of design problems that can be 
pursued with the code. For example one might wish to 
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determine the smallest length condenser section and the best 
internal geometry for a specified heat transfer rate-- 
perhaps to cool an electric motor. 

To maximize heat transfer rate through the condenser 
wall the designer adds a number of fins to increase the 
inside surface area. As more and more fins are added 
however, the cross-sectional area for conduction through 
each fin is decreased. Also, the film thickness of the 
condensate in the though increases and in fact could com- 
pletely cover the fins and substantially reduce heat trans- 
fer through the fin. So there should exist some optimum 
combination of number of fins, fin height, fin half angle, 
and ratio of trough width to fin width that will permit 
maximum heat transfer rate. 

The design study undertaken was to determine the fin 
height, fin half angle, and fin spacing which would yield 
the maximum heat transfer possible. The design variables 
then were BFIN, FANGL, and ZOA. Other potential design 
variables were held constant. The objective function to be 
maximized was QTOT , heat transfer rate out of the condenser. 

For comparison, the theoretical upper limit on heat 
transfer was calculated based on an external surface tempera- 
ture equal to the working fluid saturation temperature. This 
assumes that there is no thermal resistance across the 
condensate and the condenser wall. When this upper limit 
was used, the maximum heat transfer rate was predicted to be 
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63,322 BTU/HR using the following formula: 



Q = h 2TTr L (Twall - T°») (8) 

where 

h = outside convective heat transfer coefficient 
(5000 BTU/HR-FT 2 -°F) 

r = average outside radius of condenser wall (0.0056 ft) 
L = condenser length (0.75 feet) 

Twall=temperature of the outside wall (100°F) 

T 00 = ambient temperature (60°F) 

Certain constraints were placed on the design based on 
engineering judgement. For example, the number of fins was 
not allowed to exceed 400 and the minimum fin half angle 
allowed was 10 degrees. These values were based on structural 
and manufacturing considerations. 
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IV. RESULTS 



A. INTRODUCTION 

The computer code was used in conjunction with the COPES/ 
CONMIN program for sensitivity analysis and design optimiza- 
tion for maximum heat transfer rate. Numerical results are 
discussed below. 

B. SENSITIVITY ANALYSIS 

As shown in Figures 7 and 8, the heat transfer rate in- 
creases for an independent increase of the design variable 
of interest, but levels off at the theoretical maximum heat 
transfer rate for the given overall geometry. The heat trans- 
fer rate then is limited by the external resistance which has 
become the controlling factor. 

In Figure 9, the heat transfer rate increases in a simi- 
lar way for an increase in external heat transfer coefficient. 
Again, the rate of increase appears to lessen for external 
heat transfer coefficients above 10,000 BTU/HR-FT *°F due to 
other limiting resistances. In Figure 10 the heat transfer 
rate is observed to increase linearly with condenser length. 

This is expected since Q is a function of the area and for small 
values of <J> the area varies directly with length. In Figure 11, 
Q is seen to rise in a non-linear manner for an increase in 
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cone half angle, <p . This is because the internal heat trans- 
fer coefficient, h^ is a non-linear function of cone half 
angle. 

Purnomo [1] concluded that the heat transfer rate continu- 
ously increased as the fin half angle decreased and that this 
was largely a result of the fact that the number of fins in- 
creased at the same time. He also stated that the increase 
in heat transfer was only slight when the fin half angle was 
less than 11 degrees. This is somewhat misleading in that in 
his analysis the number of fins was being changed with every 
change in fin half angle. This author has drawn a different 
conclusion. Figure 12 shows a plot of heat transfer rate vs. 
varying fin half angle all for a condenser with 40 fins. The 
heat transfer rate as a function of fin half angle rises 
sharply from 1-11 degrees and continues to rise, but less 
steeply, as the fin half angle increases beyond 11 degrees. 

The maximum fin half angle possible with 40 fins is 67.5 
degrees . 

C. CONSTRAINED OPTIMIZATION 

In the design problem undertaken to determine the optimum 
internal geometry for maximum heat transfer, numerous runs 
were made for condensers made of copper, stainless steel, and 
a ceramic material. These materials have thermal conduc- 
tivity values of 231,9, and 1.0 BTU/HR*FT*°F respectively. 

It was expected that for each material a different optimum 
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design would emerge. The results, therefore, were unexpected. 
From different starting points (original designs) , some out- 
side the feasible region, and for external heat transfer 

2 

coefficients *from 1000-50,000 BTU/HR'FT *°F, the same optimum 
design for maximum heat transfer was reached, which is out- 
lined in Table III below. Each material, of course, had a 
different heat transfer rate even though the geometry for 
maximum heat transfer was the same. 

TABLE III 

CONSTRAINED OPTIMIZATION RESULTS FOR COPPER 
STAINLESS STEEL, AND CERAMIC MATERIAL 

Optimum Design for All Materials 



fin height 


0.023 


inches 


fin half angle 


10.0 


degrees 


number of fins 


400 . 0 




ratio of trough width 


0.5 




to fin base width 







MATERIAL 


^external 
(BTU/HR-FT 2 -°F) 


HEAT TRANSFER RATE 
( BTU/HR) 


Copper 


1000 


13,822 




50,000 


296,560 


Stainless Steel 


1000 


10,600 




50,000 


38,650 


Ceramic Material 


1000 


3844 




50,000 


5237 
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D. UNCONSTRAINED OPTIMIZATION 

The results of the constrained optimization runs were so 
unexpected that it was decided to remove the constraint on 
the number of fins and the trough width to fin width ratio 
and repeat the optimization runs. If another identical design 
was again reached for all materials, the code would have to be 
considered in error. The results are presented in Table IV. 

The optimum design for maximum heat transfer was different 
for each material and each heat transfer coefficient. For 
each material, the optimum fin height b increases as the ex- 
ternal heat transfer coefficient increases. This provides 
more fin surface area to increase heat transfer rate as the 
outside heat transfer resistance no longer controls. Despite 
a fifty-fold change in heat transfer coefficient, each 
material maintained essentially a constant 5 */b ratio. The 
higher conductivity materials required less exposed fin sur- 
face than the lower conductivity materials. Copper, for 
instance, could be 57 percent covered by the trough condensate 
while the ceramic material was only 14-19 percent covered for 
maximum heat transfer. 

In each material the number of fins decreased for an 
increase in external heat transfer coefficient. This provides 
more space in the troughs to carry the increased condensate 
which results. Otherwise, fin performance would be degraded. 
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UNCONSTRAINED OPTIMIZATION RESULTS 
FOR COPPER, STAINLESS STEEL, AND CERAMIC MATERIAL 
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igures 



For an increase in external heat transfer coefficient the 
highest percentage change in heat transfer rate occurs in the 
material with the highest thermal conductivity. In the 
ceramic material, for instance, a fifty-fold increase in ex- 
ternal heat transfer coefficient results in only a 36 percent 
increase in heat transfer rate. In the copper condenser, the 
heat transfer rate is increased by a factor of 21. This shows 
the strong dominance of the wall resistance to heat transfer 
in the ceramic material. 

For all three materials, the trough is essentially eli- 
minated for the optimum design. Why the trough is not 
entirely eliminated in the stainless steel condenser section 
is not clear. 

Because the manufacture of a very large number of fins in 
a small diameter condenser is not practical, it is desirable 
to use the constrained optimization results as a design guide- 
line . 

A comparison of the unconstrained and constrained design 
for a copper condenser section, for instance, shows that 
using the more realistic number of 400 fins and its proper 
fin height b instead of 2040 fins results in a heat transfer 
rate difference of only 0.6 percent. This is shown clearly 
in Table V. 
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TABLE V 



Comparison of Unconstrained and Constrained Optimum Design for 
a Copper Condenser Section with External^Heat Transfer Coeffi- 
cient of 50000 BTU/HR-FT -°F. 



b (inches) 

constrained 0.023 

unconstrained 0.019 



Number %Difference 

of Fins Q ( BTU/HR) in Q 

400 296,560 

2040 298,330 0.6% 



The author's conclusion is that within the realm of what 
can now be manufactured, the same basic design is best for 
all materials regardless of the external heat transfer 
coefficient. That is, the designer should machine as many 
fins as the condenser material and the manufacturing process 
will allow. The code should be used to determine the fin 
height to avoid degrading the fin efficiency by too high a 
level of condensate in the trough. As seen in Table V, it 
may be possible to lower the constraint on the number of fins 
more than once and compare the resulting heat transfer rate 
to achieve an effective but less expensive design to manu- 
facture . 

E. A CAUTIONARY NOTE 

The reader should be cautioned that when an analysis code 
based on assumptions made for certain conditions is linked to 
an optimizer, the code may change the geometry or other 
conditions such that the original assumptions are no longer 
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valid. Judicious use of side constraints and engineering 
common sense can prevent disaster. Numerical solutions must 
not be blindly accepted. 

It should be noted that all calculations in this thesis 
were made with water as the working fluid. Sensitivity data 
were run with an original rotational speed of 3000 RPM while 
design optimization was done for a rotational speed of 3600 
RPM. Different 'combinations of operating parameters, working 
fluids, and condenser section materials may not behave in a 
predictable manner. While certain trends are shown for* the 
particular condenser section studied here , it would be 
preferable to modify the code to analyze the particular 
problem at hand than to extrapolate these findings to too 
broad an application. 

While a design improvement is almost certain using an 
optimization routine like COPES/CONMIN , there is no guarantee 
that a global optimum has been found. Therefore, for a given 
design problem, the designer should start the optimization 
process from several different original designs and compare 
the results. 
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V. CONCLUSIONS 



1. The computer code is valid and can be used for single 
analysis, sensitivity analysis, and automated design of an 
internally finned rotating heat pipe. It can also be used to 
generate numerical data to find a correlation useful for 
design . 

2. For zero fins, the code converges to the results for 
a smooth tube obtained by Schafer [6]. 

3 . The heat transfer rate increases for an independent 
increase in fin half angle, rotational speed, and number of 
fins but levels off at the theoretical maximum heat transfer 
rate for the heat pipe. The heat transfer rate continues to 
increase for an increase in condenser length and cone half 
angle . 

4. Within the realm of known materials, maximum heat 
transfer occurs for the same fin geometry regardless of the 
external heat transfer coefficient. That is, for a given 
condenser radius, one should machine as many fins as possible 
with the condenser material used to maximize heat transfer. 
The code should be used to determine the correct fin height 
to avoid degradation of fin performance by too high a level 
of condensate in the trough. 
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5. When the constraint on minimum fin half angle is re- 
moved the optimum design is a large number of very thin fins. 
A practical design then might be a series of long rectangular 
fins . 
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VI. RECOMMENDATIONS 



Based upon the calculated results of this thesis, the 
following recommendations are made. 

1. Build an internally finned heat pipe to obtain 
experimental data for comparison to the analytical 
predictions . 

2. Generate a new model to analyze the case where cone 
half angle is zero. This would greatly decrease the 
cost of manufacture, and permit the use of internally 
finned heat exchanger tubing which is commercially 
available . 

3. Analyze different shaped fins, including rectangular. 

4. Obtain numerical data using the code to attempt to 
find a correlation for the heat transfer rate as a 
function of important variables: 

Q = Q (temperature, condenser geometry, fluid 

properties, condenser material properties, 
heat transfer coefficients, etc.) 

5. Modify the code to include as constraint functions 
structural failure modes such as rupture due to the 
internal spinning mass, buckling, and whirling. 
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Figure 1. 




Schematic Drawing of a Rotating Heat Pipe 



48 



Heat Transfer Rate, q BTU/hr 



1.7 



17.2 



Saturation Pressure, P g lbf/in^ 
2.9 4.7 7.5 11.5 




Figure 2. Operating Limits of a Typical, Water-Filled 
Rotating Heat Pipe 
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Figure 3. Internally Finned Condenser Geometry, 
Showing Fins, Troughs, and Lines of 
Symmetry 
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c = 0.03125 in 
b =0.025 in 
R 2 = 0.762 in 
0 = 1 ° 



Figure 4. Condenser Geometry Considered with 25 
Linear Triangular Finite Elements 
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i!l + lil = o 

3x 2 3y 2 



b • c . 

3T 

a ) - k ^ (T - T gat ) Along Boundary [1] 

t>) - k ^ (T - T a ) Along Boundary [3] 

3 T 

c) ^ = 0 Along Boundaries [2] and [4] 



Figure 5. Differential Equation and Boundary Conditions 
Considered in the Analysis of Purnomo [l] 
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SCHAFER^] 
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Figure 6 . Comparison of Trough Condensate Film 

Thickness vs. Length Along the Condenser 
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Figure 7. Heat Transfer Rate (Q) vs. Number of Fins 
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Figure 8. Heat Transfer Rate (Q) vs. RPM 
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Figure 9. Heat Transfer Rate (Q) vs. External Heat 
Transfer Coefficient 
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Figure 10. Heat Transfer Rate (Q) vs. Condenser 
Length 
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Figure 11. Heat Transfer Rate (Q) vs. Cone Half Angle 
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Figure 12. Heat Transfer Rate (Q) vs. Fin Half Angle 



APPENDIX A 



USER'S MANUAL 

This appendix describes the data cards for the use of 
the computer program. 

The data is divided into the COPES/CONMIN section and the 
heat pipe analysis program section. 

The COPES data is segmented into "blocks" for convenience. 
All formats are alphanumeric for TITLE and END cards, F10 for 
real data, and 110 for integer data. Comment cards may be 
inserted anywhere in the data deck prior to the END card and 
are identified by a dollar sign ($) in column one. The COPES 
data deck must terminate with an end card containing the word 
'END' in columns 1-3. 

Information is included in this appendix pertaining to 
data needed for single analysis, sensitivity analysis, or 
optimization using COPES/CONMIN. 

UNFORMATTED DATA INPUT 

While the user's sheet defines COPES data in formatted 
fields of ten, the data may actually be read in a simplified 
fashion by separating data by commas or one or more blanks. 

If more than one number is contained on an unformatted data 
card, a comma must appear somewhere on the card. If exponen- 
tial numbers such as 2.5+10 are read on an unformatted card, 
there must be no embedded blanks. Unformatted cards may be 
intermingled with formatted cards. Real numbers on an 
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unformatted card must have a decimal point. 
EXAMPLES 



Unformatted data; 

5, 7, 1.3, 1.0+20, -5.1 
5, 7, 1.3, 1.0+20, ,-5.1 
5 7 1.3 1.0 + 20 , , -5.1 
5 7 1.3, 1.0+20 0 -5.1 

Equivalent formatted data; 
col-*? 10 20 30 40 50 60 

5 7 1.3 1.0+20 0 -5.1 



70 80 
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DATA BLOCK A 



DESCRIPTION : COPES TITLE CARD 

FORMAT: 20A4 



1 


2 3 4 5 6 7 8 


Title 





FIELD 


CONTENTS 


1-8 


Any 80 Character title may be given on this card 
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DATA BLOCK B 



DESCRIPTION : PROGRAM CONTROL PARAMETERS 

FORMAT: 7110 



1 


2 3 4 5 6 7 8 


NCALC 


NDV NS V 



FIELD 


CONTENTS 


1 


NCALC: Calculation control 


2 


0- Read input and stop. Data o5 blocks A,B 
and V is required. Remaining data is 
optional . 

1- One cycle through program. The same as 
executing ANALI2 stand-alone. Data of 
blocks A,B and V is required, remaining 
data is optional. 

2- Optimization. Data of Blocks of A-I and V 
is required. Remaining data is optional. 

3- Sensitivity analysis. Data of blocks, 
A,B,P,Q and V is required. Remaining data 
is optional. 

NDV: Number of independent design variables in 

optimization . 


3 


NSV: Number of variables on which sensitivity 

analysis will be performed. 
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DATA BLOCK C 



Omit if NDV=0 in block B 



DESCRIPTION : INTEGER OPTIMIZATION CONTROL PARAMETERS 

FORMAT: 7110 



12345678 



IPRINT ITMAX ICNDIR NSCAL ITRM LINOBJ NACMX1 



FIELD CONTENTS 

1 IPRINT- Print control used in the optimization 

program CONMIN. 

0- No print during optimization. 

1- Print initial and final optimization 
information . 

2- Print above plus objective function value 
and design variable values at each iteration. 

3- Print above plus constraint values, direction 
vector and move parameter at each iteration. 

4- Print above plus gradient information. 

5- Print above plus each proposed design vector, 
objective function and constraint values 
during the one dimensional search. 

2 ITMAX: Maximum number of optimization inter- 

actions allowed. Default=20. 

3 ICNDIR: Conjugate direction restart parameter. 

GT . 0- Scale design variables to order of 
magnitude one every NSCAL iterations. LT . 0- 
Scale design variables according to user in- 
put scaling values. If not zero, NDV + 1 
is recommended. 

5 ITRM: Number of consecutive iterations which 

must satisfy relative or absolute convergence 
criterion before optimization process is 
terminated. Default=3. 

6 LINOBJ: Linear objective function identifier. 

If the optimization objective is known to be 
a linear function of the design variables, 
set LIN0BJ=1. Def ault=Non-linear . 
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7 



NACMX1: One plus the maximum number of active 

constraints anticipated. Def ault-NDV-2 . 
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DATA BLOCK D 



Omit if NDV = 0 in block B 



DESCRIPTION : FLOATING POINT OPTIMIZATION PROGRAM PARAMETERS 

FORMAT: 7F10 



1 


2 3 4 5 6 7 8 


FDCH 


FDCHM CT CTMIN CTL CTLMIN THETA 



NOTE: Two cards are read here. 



FIELD 


CONTENTS 


1 

2 

3 

4 

5 

6 

7 


FDCH: Relative change in design variables in 

calculating finite difference gradients. 
Default=0 . 01 . 

FDCHM: Minimum absolute step in finite dif- 

ference gradient calculations. Defaults 
0 . 001. 

CT : Constraint thickness parameter. Defaults 

=-0.05. 

CTMIN: Minimum absolute value of CT consi- 

dered in the optimization process. 
Default=0 . 004 . 

CTL: Constraint thickness parameter for 

linear constraints. Default=-0 . 01 . 

CTLMIN: Minimum absolute value of CTL con- 

sidered in the optimization process. 
Default = 0 . 001 . 

THETA: Mean value of push-off factor in the 

method of feasible directions. Defaults 

1.0 . 
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DATA BLOCK D 



Omit if NDV=0 in block B 



FORMAT : 


4F10 


1 


2 3 4 5 6 7 8 



DELFUN DABFUN ALPHAX ABOBJ1 



FIELD 


CONTENTS 


1 

2 

3 

4 


DELFUN: Minimum relative change in objective 

function to indicate convergence of the 
optimization process. DEFAULT=0 . 001 . 

DABFUN: Minimum absolute change in objective 

function to indicate convergence of the 
optimization process. DEFAULT=0 . 001 times 
the initial objective value. 

ALPHAX: Maximum fractional change in any design 

variable for first estimate of the step in 
the one-dimensional search. DEFAULT=0.1. 

AB0BJ1: Expected fractional change in the 

objective function for first estimate of the 
step in the one-dimensional search. 

DE FAULT =0 . 1 . 


REMARKS 





1) The DEFAULT values for these parameters usually work well. 
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DATA BLOCK E Omit if NDV=0 in block B 



DESCRIPTION: 


TOTAL NUMBER OF DESIGN VARIABLES, DESIGN 
OBJECTIVE IDENTIFICATION AND SIGN 



FORMAT : 2110, F10 



1 


2 3 4 5 6 7 8 


NDVTOT 


IOBJ SGNOPT 



FIELD 


CONTENTS 


1 


NDVTOT: Total number of variables linked to the 


2 


design variables. This option allows two or 
more parameters to be assigned to a single 
design variable. The value of each para- 
meter is the value of the design variable 
times a multiplier, which may be different 
for each parameter. DEFAULT=NDV. 

IOBJ: Global variable location associated with 


3 


the objective function in optimization 
SGNOPT: Sign used to identify whether function 

is to be maximized or minimized. +1.0 indi- 
cates maximization. -1.0 indicates minimi- 
zation. If SGNOPT is not unity in magnitude, 
it acts as a multiplier as well, to scale the 
magnitude of the objective. 
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DATA BLOCK F Omit if NDV=0 in block B 



DESCRIPTION: 


DESIGN VARIABLE BOUNDS, INITIAL VALUES AND 
SCALING FACTORS 



FORMAT: 4F10 



1 


2 3 4 5 6 7 8 


VLB 


VUB X SCAL 



NOTE: READ ONE CARD FOR EACH OF THE NDV INDEPENDENT DESIGN 

VARIABLES. 



FIELD 


CONTENTS 


1 

2 

3 

4 


VLB: Lower bound on the design variable. If 

VLB . LT . -1 .OE+15 , no lower bound. 

VUB: Upper bound on the design variable. If 

VUB . GT . 10 .E+15 , no upper bound. 

X: Initial value of the design variable. If 

X is non-zero, this will supercede the value 
initialized by the user-supplied subroutine 
ANALIZ. 

SCAL: Design variable scale factor. Not used 

if NSCAL.GE.O in BLOCK C. 
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DATA BLOCK G Omit if NDV=0 in block B 
DESCRIPTION : DESIGN VARIABLE IDENTIFICATION 

FORMAT : 2110, F10 



1 


• 

2 3 4 5 6 7 8 


NDSGN 


IDSGN AMULT 



NOTE: READ ONE CARD FOR EACH OF THE NDVTOT DESIGN VARIABLES 

IN THE SAME ORDER AS IN BLOCK F. 



FIELD 


CONTENTS 


1 

2 

3 


NDSGN: Design variable number associated with 

this variable. 

IDSGN: Global variable number associated with 

this variable. 

AMULT: Constant multiplier on this variable. 

The value of the variable will be the value 
of the design variable, NDSGN, times AMULT. 
DEFAULT = 1. 0 . 
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DATA BLOCK H Omit IF=NDV 0 in block B 
DESCRIPTION : NUMBER OF CONSTRAINED PARAMETERS 

FORMAT: 110 



1 


2 3 4 5 6 7 8 


NCONS 





FIELD 


CONTENTS 


1 


NCONS: Number of constraint sets in the 

optimization problem. 


REMARKS : 





1) If two or more adjacent parameters in the global common 
block have the same limits imposed, these are part of the 
same constraint set. 
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DATA BLOCK: I 



Omit if NDV=0 in block B, or NCONS=0 in 
block H 



DESCRIPTION: 


CONSTRAINT IDENTIFICATION AND CONSTRAINT BOUNDS 


FORMAT: 3110 




1 


2 3 4 5 6 7 8 



ICON JCON LCON 

NOTE: READ TWO CARDS FOR EACH OF THE NCONS CONSTRAINT SETS. 



FIELD 


CONTENTS 


1 

2 

3 


ICON: First global number corresponding to the 

constraint set. 

JCON: Last global number corresponding to the 

constraint set. DEFAULT=ICON. 

LCON: Linear constraint identifier for this 

constraint set. LCON=l indicates linear 
constraints . 
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DATA BLOCK: I 



(CONTINUED) 



FORMAT : 4F10 

12345678 



BL SCAL1 BU SCAL2 



FIELD CONTENTS 



1 

2 

3 

4 



BL: Lower bound on the constrained variables. 

If BL.LT.-1.0E+15 , no lower bound. 

SCAL1: Normalization factor on lower bound. 

DEFAULT=MAX of ABS(BL), 0.1. 

BU: Upper bound on the constrained variables. 

If BU. GT . 1 . OE+15 , no upper bound. 

SCAL2 : Normalization factor on upper bound. 

DEFAULT =MAX of ABS(BU), 0.1. 



REMARKS : 

1) The normalization factor should usually be defaulted. 

2) The constraint functions sent to CONMIN are of the form: 

(BL- VALUE )/SCALl . LE . 0.0 and (VALUE - BU)/SCAL2 . LE . 

0 . 0 . 

3) Each constrained parameter is converted to two con- 

constraints in CONMIN unless ABS(BL) or ABS(BU) exceeds 
1. 0E + 15, in which case no constraint is created for that 
bound . 
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DATA BLOCK : P Omit if NSV=0 in block B 

DESCRIPTION : SENSITIVITY OBJECTIVES 

FORMAT: 2110 



12345678 



NSOBJ IPSENS 



NOTE: TWO OR MORE CARDS ARE READ HERE. 



FIELD CONTENTS 

1 NSOBJ: Number of separate objective functions 

to be calculated as function of the sensi- 
tivity variables. 

2 IPSENS: Print control. If IPSENS. GT.O, de- 

tailed print will be called at each step in 
the sensitivity analysis. DEFAULT=No 
print . 
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DATA BLOCK: P 



(CONTINUED) 



DESCRIPTION : 

FORMAT : 8110 

1 2 3 4 5 6 7 ' 8 



NSN1 NSN2 NSN3 NSN4 



FIELD CONTENTS 

1-8 NSNI: Global variable number associated with 

the sensitivity objective functions. 

REMARKS : 

1) More than eight sensitivity objectives are allowed. Add 
data cards as required to contain data. 
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DATA BLOCK : Q_ Omit if NSV = 0 in block B 

DESCRIPTION: SENSITIVITY VARIABLES 

" / 

FORMAT: 2110 



12345678 



ISENS NSENS 



NOTE: READ ONE SET OF DATA FOR EACH OF THE NSV SENSITIVITY 

VARIABLES. TWO OR MORE CARDS ARE READ FOR EACH SET 
OF DATA. 



FIELD CONTENTS 

1 ISENS: Global variable number associated with 

the sensitivity variable. 

2 NSENS: Number of values of this sensitivity 

variable to be read on the next card. 
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DATA BLOCK: Q 



(CONTINUED) 



DESCRIPTION: 
FORMAT : 


8F10 


1 


2 3 4 5 6 7 8 



SNS1 SNS2 SNS 3 SNS4 



FIELD 


CONTENTS 


1-8 


SENSI : Values of the sensitivity variable . 

I=1,NSENS. 1=1 corresponds to the nominal 
value . 


REMARKS : 





1) More than eight values of the sensitivity variable are 

allowed. Add data cards as required to contain the data. 
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DATA BLOCK: V 



DESCRIPTION : COPES DATA 'END' CARD 

FORMAT : 3A1 

1 2 3 4 5 6 7 8 



END 



FIELD CONTENTS 

1 The word 'END' in columns 1-3. 

REMARKS 

1) This card MUST appear at the end of the COPES data. 

2) This ends the COPES input data. 

3) Data for the user-supplied routine, ANALIZ, follows this. 
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HEAT PIPE ANALYSIS 



Data for the heat pipe analysis follows the 'END' card 
in the COPES data deck. If - the general design capability of 
COPES/CONMIN is not needed, the heat pipe analysis can 
be run by setting NCALC = 1 in field number 1 of data block 
B; or in a stand-alone mode by using the following main 
program. 



C MAIN PROGRAM FOR HEAT PIPE ANALYSIS 
C READ, EXECUTE, AND PRINT RESULTS 
DO 10 ICALC = 1,3 
10 CALL ANALIZ (ICALC) 

STOP 

END 



* 



79 



DATA BLOCK: AA 



DESCRIPTION: 
FORMAT: 315 


ELEMENT CONNECTIVITY 


1 


2 3 4 5 6 7 8 


NEL 


NSNP NBAN 



FIELD 


CONTENTS 


1 

2 

3 


NEL: Number of elements 

NSNP: Number of system nodal points 

NBAN : System band width 
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DATA BLOCK: BB 



DESCRIPTION : ELEMENT CONNECTIVITY 

FORMAT: 415 I = 1,3; IEL=1, NEL 



1 


2 3 4 5 6 7 8 


IEL 


IC0R( IEL , I ) 



FIELD 


CONTENTS 


1 

2 


IEL: The element number 

ICOR(IELjI): System nodal point corresponing 

to nodal point I of element IEL 


REMARKS : 





1) Number all elements with convective boundaries first from 
top to bottom, then number the remaining elements. 

2) Number the nodal points of each element moving in a 
counter-clockwise direction. 

3) The elements with convective boundaries have nodal points 
1 and 2 located on the convective boundary . 
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DATA BLOCK: CC 



DESCRIPTION : CONDENSER GEOMETRY 

FORMAT: 7G10 . 5 



1 


2 3 4 5 6 7 8 


CLI 


CANGL RBASEI R2I THICKI BFIN TZ 



FIELD 


CONTENTS 


1 

2 

3 

4 

5 

6 
7 


CLI: Condenser length (inches) 

CANGL: Cone half angle (degrees) 

RBASEI: Inside radius of condenser base 

( inches ) 

R2I: Intermediate radius (inches) 

THICKI: Condenser wall thickness (inches) 

BFIN: Height of fin (inches) 

* TZ : Nodal point temperature initial guess 

(degrees F) 


REMARKS 





] ) Set TZ equal to TSS-Saturation temperature of the 
working fluid 
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DATA BLOCK: DD 



DESCRIPTION : FINITE ELEMENT GEOMETRY 

FORMAT : 515 

12345678 



NDIV NEST NEFB NBOTI NBOTF 



FIELD CONTENTS 

1 NDIV: Number of increments along the length 

of the condenser. 

2 NEST: Number of the element on the right end 

of the trough. 

3 NEF3: The element number with convective 

boundary located at the base of the fin. 

4 NBOTI: The element number with convective 

boundary located at the right hand of the 
bottom side. 

5 NBOTF: The element number with convective 

boundary located at the left hand of the 
bottom side. 
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DATA BLOCK: EE 



DESCRIPTION: 


DATA FOR RUNNING 


FORMAT: 4F10 


.2 


1 


2 3 4 5 6 7 8 


RPM 


TSS TINF HINF 



FIELD 


CONTENTS 


1 

2 

3 

4 


RPM: Rotation rate of heat pipe (RPM) 

TSS: Saturation temperature of the working 

fluid (degrees F) 

TINF: Outside temperature (degrees F) 

HINF: Outside convective heat transfer 

coefficient (3TU/HR-FT /o F) 
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DATA BLOCK: FF 



DESCRIPTION: 


CONVERGENCE CRITERION 


FORMAT: G10 . 


.9 


1 


2 3 4 5 6 7 8 


CRIT 





FIELD 


CONTENTS 


1 


CRIT: Convergence criterion on finite element 

solution for incremental heat transfer rate 
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DATA BLOCK: GG 



DESCRIPTION: 


INTERNAL FIN GEOMETRY 


FORMAT: 2G10 


.5 


1 


2 3 4 5 6 7 8 


FANGL 


ZOA 



FIELD 


CONTENTS 


1 

2 


FANGL: Fin half angle (degrees) 

ZOA: Ratio of trough width to fin base width 
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DATA BLOCK: HH 



DESCRIPTION: 
FORMAT: 15 


INTERNAL FIN GEOMETRY 


1 


* 2 3 4 5 6 7 8 


IFF 





FIELD 


CONTENTS 


1 


IFF: (n-1), where n is the number of rows of 

the upper triangular fin section 



NOTE: See Figure A'-l. 
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DATA BLOCK: II 



DESCRIPTION : INTERNAL FIN GEOMETRY 

FORMAT : 1615 1=1, IFF 

12345678 



KFIN(I) KFF(I) 



FIELD CONTENTS 

1 KFIN : The number of system nodal points located 

on the symmetric boundary of triangular fin 
section, but does not include the system 
nodal points located at the base of the fin 
and the apex. 

2 KFF : The number of system nodal points located 

along the fin convective boundary, but does 
not include the system nodal points located 
at the base of the fin and the apex. 



NOTE: See Figure A-l. 
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DATA BLOCK: JJ 



DESCRIPTION : INTERNAL FIN GEOMETRY 

FORMAT : 2G10.5,4I5 



1 


t 

2 3 4 5 6 7 8 


DOBF 


DOTH JTC JLC JINT KT 



FIELD 


CONTENTS 


1 

2 

3 

4 

5 

6 


DOBF: Number of column within fin. 

DOTH: Number of column within though. 

JTC: The number of the system nodal point 

located at the junction of the symmetry 
boundary and the line of intersection 
between the fin and the condenser wall. 

JLC: The number of the system nodal point 

located at the center of system coordinates . 

JINT: The numerical difference between the two 

adjacent system nodal points vertically 
at the condenser section. 

KT: The number of rows within the wall section 



NOTE: See Figure A-l . 
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IFF 


= 2 




DOBF 


KFIN 


= 2 ; 


4 


DOTH 


KFF 


= 3 ; 


6 


JTC 



= 3.0 


JLC = 17 


= 1.0 


JIN'? = 5 


= 7 


KT = 2 



Figure A-l. Specification for Input Data to Determine 
the Coordinates of the System Nodal Points 
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noonoooooooooooooono 



APPENDIX B 



COMPUTER PROGRAM LISTING 



***************$$$*******#$***#*£*«****************;):**$* 
* 

* ANALYSIS CF ROTATING HEAT PIPE , USING TRIANGULAR 

* ELEMENT MODEL 

* COMPILED BY MAJOR IGNATIUS.S. PUR NO MO IN JUNE 1978 

* 

* MODIFIED TO PERMIT NUMERICAL OPTIMIZATION 

* USING COPES/CONMIN 

* BY LCDR WILLIAM A. DAVIS, JR. IN SEPTEMBER 1980 

* 

$#****$*#$$***:*******:*$*$:****:*** ********* * ****#***$:(:** $* 



SUBROUTINE ANALIZ CHANGES THE ARRAY OF DESIGN 
VARIABLES FROM SINGLE TO DOUBLE PRECISION AND 
BACK. COPES/CONMIN USES SINGLE PRECISION ONLY; 
DOUBLE PRECISION IS MAINTAINED IN SUBROUTINE FUN 
TO ALLOW FOR POSSIBLE ILL CONDITIONING. 



10 

20 

30 



40 



C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUBROUTINE ANALIZ (ICALC) 

COMMON / GLCBCM/ ARRAY (750) 
COMMON /GLCB1/ BARAY (50) 

REALMS BARAY 

IF ( ICALC. GT.l) GO TO 20 

DO 10 1=1,50 

BARAY ( I ) = 0 .0 DO 

CONTINUE 

CONTINUE 

CO 30 1=1,50 

BARA Y ( I ) =DBLE (ARRAY ( I ) ) 

CALL FUN (ICALC) 

DO 40 1=1,50 

ARRA Y( I ) = SNGL (BARA Y ( I ) ) 

CONT INUE 

RETURN 

END 

GUIDE TO FORTRAN VARIABLE NAMES 



ALFA FIN HALF ANGLE (RADIANS) 

6F I N HEIGHT CF FIN (INCHES) 

B V IN HEIGHT OF FIN (FEET) 

CALFA COSINE OF ALFA 

CANGL CONE HALF ANGLE (DEGREES) 

CBASE INSIDE CIRCUMFERENCE OF CONDENSER (FEET) 

CEXIT INSIDE CIRCUMFERENCE AT CONDENSER EXIT (FEET) 

CF THERMAL CONDUCTIVITY OF CONDENSATE FILM (BTU/HR) 

CL CONDENSER LENGTH (FEET) 

CLI CCNCENS ER LENGTH (INCHES) 

C PHI COSINE OF PHI 

CRIT CONVERGENCE CRITERION 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

£ 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



OIV 

DMTOT 

DOBF 

DOTH 

FANGL 

H 

HFG 
I EL 
JLC 

JTC 



KFF 
KF IN 



KT 

NB AN 

NBOTF 

NBOTI 

NDIV 

NEFB 

NEL 

NEST 

NSNP 

PHI 

PI 

RBASE 

RBASEI 

REXIT 

SALFA 

SPHI 

THICK 

THICKI 

TPHI 

ZF IN 

ZOA 



FLOATING POINT VALUE OF NDIV 
CONDENSATE MASS FLOW RATE 
NUMBER OF COLUMN WIWTHIN THE FIN 
NUMBER OF COLUMN WITHIN THE TROUGH 
FIN HALF ANGLE (DEGREES) 

CONVECTIVE HEAT TRANSFER COEFFICIENT ( BTU/HR FT 2) 
LATENT HEAT OF VAPORIZATION ( BTU/L8M) 

THE ELEMENT NUMBER 

NUMBER OF SYSTEM NODAL POINT LOCATED AT 
THE CENTER OF SYSTEM COORDINATE 
NUMBER OF SYSTEM NODAL POINT LOCATED AT 
THE JUNCTION OF THE SYNNETRY BOUNDARY AND 
THE LINE OF INTERSECTION BETWEEN THE FIN 
AND THE CONDENSER WALL 

NUMBER OF SYSTEM NODAL POINTS LOCATED ALONG 
THE FIN CONVECTIVE BOUNDARY 
NUMBER OF SYSTEM NODAL POINTS LOCATED CN THE 
SYSMMETRIC BOUNDARY OF TRIANGULAR FIN SECTION 
NOT COUNTING POINTS AT BASE AND APEX 
NUMBER OF ROWS WITHIN THE WALL SECTION 
SYSTEM BAND WIDTH 
LAST ELEMEMT AT BOTTCM SIDE 
FIRST ELEMENT AT BCTTGM SIDE 
NUMBER OF INCREMENT 
ELEMENT NUMBER AT BASE OF FIN 
NUMBER OF ELEMENTS 
ELEMENT NUMBER AT END OF TROUGH 
NUMBER OF SYSTEM NODAL POINTS 
CONE HALF ANGLE (RADIANS) 

PI 

INSIDE RADIUS OF CONDENSER BASE (FEET) 

INSIDE RADIUS OF CONDENSER BASE (INCHES) 

INSIDE RADIUS OF CONDENSER EXIT (FEET) 

SINE OF ALFA 
SINE CF PHI 

CONDENSER WALL THICKNESS (FEET) 

CONDENSER WALL THICKNESS (INCHES) 

TANGENT OF PHI 
NUMBER OF FINS 

RATIO CF TROUGH WIDTH TO FIN BASE WIDTH 



SUBROUTINE FUN READS INPUT DATA, PERFORMS HEAT 
TRANSFER ANALYSIS, AND PRINTS RESULTS. 
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SUBROUTINE FUN (ICALC) 

IMPLICIT REAL*8(A-H,0-Z) 

COMMON / GLCB1/ CL I , CANGL, RBA SE I , R2 1 , TH I CKI , BF IN, TZ , TSS 
1,TINF,HINF,FANGL,Z0A,ZFIN,B0A,QT0T,RPM 
01 MENS ION Z( 200), EPS ( 200) , HZ ( 200) ,XC0F< 5) ,COF (5) , T(200 
1 ),QB(200) tOMDCT (200 ) ,UF(200) , CF( 200 ) , C W ( 200 ) ,AMT0T ( 200 
2)»R(200) » Cl NC ( 200 ) » TB ( 200 ) »TT ( 200 ) , TI B (200 ) , QT I NC ( 200 ) 
3»NC(200) , Q TOTAL ( 100 ) » TE( 200) ,R00 TR ( A ) , ROOTI (A) , DEL (200 
4) ,TBM<200) » RHGF (200 ) 

COMMON /APOL/ OOBF , DOTH , KF IN ( 50) , KFF ( 50) , I FF , JTC, J LC, J 
1 INT, KT 

COMMON /MAFO/ A ( 200 , 50 ) , F ( 200 , 1 ) , H (200 ) , TS ( 200 ) , TS AT , C 
lK»NEL»NSNP»NBANt ICOR( 200 t 3) 

COMMCN / PCRD/ X ( 200 ) ? V ( 200 ) , EZERO, BV I M , TH ICK , TALF A , APS 
HDENIAl,Bl»ZZ)*(-1.0D0*(Al*ZZ**3/3.000+Bl*ZZ**2/2 .000 ) 

L) 

IF I CALC=1 READ INPUT CATA 
IF (ICALC. GT. I) GO TO 10 






INPUT MODE 






ELEMENT CONNECTIVITIES 

READ ( 5» 420) NEL ,N SNP , NBAN 
WRITE ( 6 1 430 ) N EL » NSNP , NBAN 

READ (5,440) ( I EL , ( I COR ( I EL , I ) , I = L , 3 ) , I EL=1 , N EL ) 

WRITE <6,450) 

WRITE (6,251) (I EL, < I COR (I EL, I) , 1=1,3) , IEL = 1,NEL) 

THE CONDENSER GEOMETRY 

READ <5,460 CLI , CANGL ,RBASEI ,R2 I ,THI C KI , BF I N,TZ 
WRITE <6,470) CL I , CANGL, RB ASE I , R 21 , TH I CK I , BF I N , TZ 
READ <5,480) NO IV, NEST ,NEFB» NBOT I , NBOT F 
WRITE <6,490) NDI V, NE ST, NE FB , NBO TI ,NBO TF 

DATA FOR RUNNING 

READ <5,500) RPM, TS S , T INF , HINF 
WRITE <6,510) RPM, TSS, TINF , HINF 

THE CONVERGENCE CRITERIAN 

READ <5,520) CRIT 
WRITE <6,530) CRIT 

INTERNAL FIN GEOMETRY 

READ <5,540) FANGL , ZO A 
WRITE <6,550) FANGL, ZOA 
READ < 5, 560 IFF 
WRITE <6,570) IFF 

READ ( 5,580) ( KF I N < I ) , KFF < I ) , I =1 , 1 FF ) 

READ (5, 590) DOBF, DO TH , J TC , JLC , J INT, KT 
NHB=NEFB/2 

NTM=NBOTI+ ( NBOTF-NBOTI ) /2 
NBF=NSOT F + l 

WRITE <6,600) IC OR (NBOT I ,2) »I COR (NEFB»1),IC0R<NTM,2), I 
1CQR(NEST, 1) , ICOR( NB OTF , 1 ) 

RETURN 

IF I CA LC = 2 PERFORM THE HEAT TRANSFER ANALYSIS 
10 IF (ICALC. GT. 2) GO TO 360 






EXECUTION MODE 
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CONVERT UMTS OF ALL DIMENSIONAL PARAMETERS 
TO FEET. CONVERT UNITS OF ANGLES TO RADIANS. 

CL=CLI / 12 .000 
R2=R 21 /12.CD0 
RBAS E=RBAS E 1/ 12.0D0 
BVIN=BFIN/12.0D0 
01 V=DFL0AT { NO I V ) 

P 1=3. 141 5926535 8979 DO 

PH I =2. 0D0*CANGL*PI /360.0D0 

SPHI = DS INI PHI ) 

CPHI =DCCS ( PH I ) 

TPHI =DTAN ( PH I ) 

CELX=CL/CIV 
CBASE=2.0D0*PI*R8AS6 
REXIT=RBASE+CL*TPHI 
CEXIT=2.000*P I*REXIT 
THICK=THICKI/12.0D0 
ALFA=FANGL*2.0D0*PI/360. 0D0 
SALFA=DS I M ALFA) 

CALFA=OCGS ( ALFA) 

TALFA=DTAN< ALFA) 

EZER0=2. ODO* 8VIN*T ALFA 

80UNCARY CCN D IT IONS AND TEMPERATURE ESTIMATES 

ALONG THE FIN BOUNDARY 

DO 20 NTINF«NBOTI, NBOTF 
20 TS (NTI NF ) =TI NF 
DO 30 NNT=NBF,NEL 
TS (NNT ) =0 . ODO 
30 H( NNT) =0. ODO 

DO 40 IGT-l.NEST 
IE=ICOR( IGT ,2 ) 

40 T( I E )= TZ 

IG=I COR ( NEST » 1 ) 

T ( IG)=TZ 

CMEGA=RPM*2.0D0*P 1*60 . ODO 
DO 50 KL=NBCTI, NBOTF 
50 H(KL)=HINF 
HI FN=H INF 
TSAT=TSS 
EPSO=ZOA*EZERO 
EOA=BV I N/ ( EZ ERO/2 . ODO) 

ZFIN=C8ASE/ (EZERO+EPSG) 

SURFAR=ZFIN*(2.0D0*( B V IN /CALFA ) +EP SO ) 
EPSEX=(CEXIT-<ZFIN*EZERO) )/ZFIN 
BETA=( EPSEX-EPSO) /DI V 
ZZERO= BVIN/CALFA 
ZA =0. ODO 

DO 60 NSAT=1,NEST 
60 T S ( NSAT ) =TSAT 

TSOLID=( TSAT+TINF) /2.0D0 
TEMPORARY CHANGE - TFILM 
CT=0. ODO 
QBTO T= 0. CDC 
QT 1=0. ODO 
QTF=C. ODO 
QT0T=0 .ODO 
DMT0T=0 .ODO 
NK=NDI V+l 
DO 350 NI=1,NK 
R(NI) = R2 + M*DELX*SPHI 
EPSINI ) = E P SC +N I *B E T A 
APS=EPS (N I ) 
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NCCAL POINT COORDINATES 



70 



CALL COORO 
Z( 1 ) =ZA 

00 70 I ZEL=1 » NEFB 
NA= I COR ( I Z EL » I) 

NB=I COR ( I ZEL 1 2 ) 

XE=X(NA)-X(NB) 

YE=Y (N A) -Y (NB) 
ELZ=DSQRT(XE**2+YE**2) 

Z( IZEL+1J=Z( IZED + ELZ 

XZB=X ( ICCPINHBt 1) )-X ( IC0R( 1, 2) ) 
YZB=Y( ICCR(NHBtl) ) -Y ( I COR ( It 2 ) ) 
ZB=DSQRT <XZB**2+YZB**2) 

ZC=ZZER0 

1 M= 1 



PARABOLIC TEMPERATURE DISTRIBUTION ALONG THE FIN 
BOUNOARY, USING LAGRANGE INTERPOLATION 



) 

) 



) +AY*TS(NY) )/ (2.0D0*AY) 

PROPERTIES 
HFG=1097.200— 0.601875DC*TS(1) 

RHOF (N I >=62.77400-0 .00 25569800 *TF-0. 00005 3572D0*TF**2 
CFINI ) =0.3 03400+0. 000 738 92 7D0*TF-0. 00 000 1473 2 100*TF**2 
UF(N 11=0.0013970 0-0. 0000 1466 9D0*TF + 0. 00000006 31253Q0*T 
1 F* *2-0 .00000 000009 7656900*TF**3 
LF(NI) =36C0*UF(NI ) 

CW (NI)=231.777200— 0.0222 20 0*T SOL ID 
CW(NI) =8. 776+0. 0026 500 *TSOLID 
CWINI >=1.0 
CW (N I J = 20000 .0 
CK=C W( NI ) 

CO NS T= RHOF (NI)**2*0MEGA**2*HFG*CPHI*CALFA*R(NI) 

INITIAL FILM THICKNESS 
0EL< 11=0.0000675200 
IF (NI .GT.l) GO TO 100 

0EL( 1) = 1. 1C7*< (.( TSAT-TINF) *CF (NI ) / ( UF ( NI ) *HFG*3600 . 000 
1 ))**0.25)*( (UF(NI)/(RHOF(NI)*OMEGA) )**0.5) 

100 CONTINUE 

AVERAGE ELEMENT CONVECTIVE COEFFICIENT ALONG 
THE FIN BOUNOARY 

ZSTAR=ZZERC-DEL (NI J/CALFA 
AZZ=DEL( NI ) /SALFA 
ZZ= ZSTAR 

AZ S=DA8S( 4*CF (NI)*UF(NI)*H0EN(A1, 81, ZZ) /CONST ) **0 .2500 
HAC=0. 000 



80 TP 1=T ( ICCR ( 1 1 2) ) 
TP2=T ( I COR (NHB» 1) ) 
TP3=T( ICCR(NEFB,1 ) ) 
AP 1= TP 1/ ( ZB*ZC ) 
AP2=TP2/ (ZB*( ZB-ZC ) 
AP3=TP3/(ZC*( ZC-ZB) 
BP 1=— ( ZB+ZC ) *AP 1 
BP2=-ZC*AP2 
BP 3=-Z8*AP3 
A1=AP1+AP2+AP3 
B1=BP1+8F2+BP3 
TC=0. 000 
CO 90 NY=ltNEST 
90 TC=TC+T ( ICCR ( NY »2 ) ) 
AY=DFLOA T( NY+ 1 ) 

TF=( TC+T ( ICCR (NY » 1) 



SOLID-FLUIO 
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00 190 I EL = 1 1 NEFB 
AZ=Z( I EL) 



BZ=Z(I EL+li 

IF (ZSTAR. LE.BZ) GO TO 110 



GO TO 120 

110 IF (HAC.NE. 0.000) GO TO 180 
BZ=ZSTAR 

120 IF (IEL.NE.1) GO TO 130 
AK=(BZ-AZ) /5.0D0 
ZZ=AK 
GO TO 140 

130 AK=(8Z-AZ) /4.0D0 
2Z=AZ 

140 ZEL=4*AK 

00 150 Nh= 1 » 5 

HZ(NH)=DABS(CF(NI)**3*C0NST/(4*UF(NI ) *HOEN ( A 1 , B 1, ZZ ) ) ) 
1**0.2500 
150 ZZ = ZZ+AK 

CGNH=AK*(HZ( 1)+4*HZ< 2) +2*HZ( 3) +4*HZ ( 4 ) +HZ ( 5 ) >/( 3*ZEL) 
IF (ZSTAR. EQ.BZ) GO TO 160 
H( IEL)=CCNH 
GO TO 190 
160 AZ=ZSTAR 

HAZ=CCNH*(AZ-Z( I EL ) > 

OELA =AZS 
170 BZ=Z ( I EL+ 1 ) 

DELB=( 8Z- ZSTAR) *AZZ/ ( ZZERO-ZSTAR ) 
0ELZ=(DELA+D£LB>/2.000 
FAC= (BZ-AZ)*CF(N I )/DELZ 
H( IEL) = ( HAZ+H AC ) / ( 8Z-Z ( IEL ) ) 

GO TO 19C 
180 AZ=Z ( I EL ) 

OELA=OELB 
HAZ=0. 000 
GO TO 170 
190 CONTINUE 

NET I=N EFB + 1 
DO 200 IEL=NETI t NEST 
200 H( IEL)=CF(NI)/DEL(NI ) 



C 

C 

c 



ENTRY INTO THE FINITE ELEMENT SOLUTION 



CALL FCRMAF 

CALL BANOEC ( NSNP, NB AN , l ) 



C 

C 

c 



THE TEMPERATURE DISTRIBUTION 




CO 210 NT= 1 , NSNP 
T ( NT)=F( NT ,1 ) 

TIB(NI ) = T( ICOR ( NBOT I ♦ 2) ) 



DO 220 NS= 1 « NSNP 




C 

C 

C 



Q AT THE BOTTOM SIDE 
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230 



240 

250 

260 



270 

280 

290 

300 

310 



320 



330 



YB=Y(NKA)-Y(NKB) 

ELB=DSQRT(XB**2+YB**2 ) 

GBI=QBI+(T(NKA)+T(NKB)-2*TS( IBEL) ) *ELB *H l IBED/2.0 00 
CB(NI)=QBI*DELX 

ITERATION UNTIL CONVERGENCE CRITERIA IS MET 

IF < IM.EQ.l) GO TO 240 

QJ=QB I 

GO TO 250 

QI=QBI 

IM=2 

GO TO 80 

AQ=DABS( CJ-CD/QJ 
IF (AQ.LE.CRIT) GO TO 260 
Cl =QJ 
GO TO 80 

OMDOTINI )=2.*QBI*DELX/HFG 
CMTQ T=0MTCT+DMD0T ( N I ) 

C1=RHCF(M )**2*OMEGA**2*R(NI )*SPHI/(3*UF(NI ) ) 

XCOF ( 1 ) =-0FT0T 
XCOF(2)=O.CCO 
XCOF ( 3 ) =0 . 000 
XC0F< 4) = C1*EPS(NI ) 

XCOF (5 ) = C1*T ALFA 



M=4 

CALL DPOLRT <XC0F,C0F, 
IF (RGOTR(l) .GT .0. 0 DO ) 
IF <R00TR<2) .GT. 0.000) 
IF (RG0TR<3) .GT.0.0D0) 
IF (RC0TR(4) .GT. 0.000) 
WRITE (6,610) 

WRITE (6,620) <R00TR< I) 
GO TC 6100 



MtROCTR, ROOTI ,IER) 
GO TO 270 



GO 

GO 

GO 



TO 

TO 

TO 



280 

290 

300 



1= 1,4) 



THE CGNOENS AT E THICKNESS 



SS4 



NI + 1 )=ROOTR( 1 ) 
G 310 
DEL(NI+1)=R00TR(2) 



GO 



GO 

DEL 



0 310 



OEUNI+l J=RCGTR(3) 



0 310 

NI+1 }=R00TR(4) 
0EL=0. 000 
IF (NI.NE.l) 



GO TO 320 



Q FROM THE TCP SIDE 

Q THROUGH FIN 

DO 330 I CEl = l ,3 
KA=IC0R(ICEL,1) 

KB=ICOR( ICEL, 2) 

XQEL=X(KA)-X(KB) 

YQEL=Y ( KB)-Y( KA ) 

ELM=DSQRT (XQEL**2+YQ EL**2) 
QEL=QEL+ (2*TS(IQEL)-T (KA)-T(KB) 
CONTINUE 

CINC(NI)=QEL*DELX 
AMTOT ( NI ) =DMTOT 
CET=GEL*CELX*ZFIN*2 
QT=QT+QET 

QA=QBI*DELX*ZFIN*2 

QTOT=QTOT+QA 



)*ELM*H( IQEL )/2.000 



C THROUGH TROUGH 
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DO 340 I QEL = 4» 4 
KA=IC0R(ICEL»1) 

KB=ICGR(IQEL»2) 

XQEL=X (KA)-X(KB) 

YQEL=Y (K8)-Y (KA) 

ELM=DSQRT ( XGEL**2+YQE L**2 ) 

CTRF=(2*TS ( IQEL )-T (KA)-T(KB) ) *ELM*H( I 3 EL ) /2. 000 
340 CONTINUE 

CTINCt NI )=QTRF*DELX 
CTCTAL(NI)=CINC(NI ) +QTINC(NI ) 
QTRFT=QTRF*DELX*ZFI N*2. 

CTF=QTF+GTRFT 
350 CONTINUE 
RETURN 

360 CONTINUE 



370 



GT I NC( NR) , QTOTAL (NR ) 



SURF AS 



380 



390 



400 



410 



420 

430 



***** OUTPUT MODE ***** 

WRITE (6,630) 

CO 370 NR= 1, NDIV 
WRITE (6,640) NR, Q I NC ( NR ) 

IwRITE ( 6,650) QT,QTF 
WRITE (6,660 ) CL I, CANGL, RBASE I.R2I , THICK I , BF I N ,RPM , TSS 
1 ,TINF,HINF,CRIT,FANGL,ZOA, IFF 
WRITE (6,670) BOA , ZOA , ZFI N ,B VI N , 

WRITE (6,680) 

00 380 NP = 1 , NSNP 

WRITE (6,690) NP , X ( NP ) , Y( NP ) , T( NP ) 

WRITE (6,700) 

CO 390 KKL = 1,N80TF 
NKX= ICOR (KKL , 1 ) 

NKY=IC0R(KKL,2) 

XP=X(NKX)-X( NKY ) 

YP=Y (NKX)-Y (NKY) 

EXY=OSQRT(XP**2+YP**2) 

CEP=DABS ( (T(NKX)+T(NKY)-2*TS(KKL ) ) *EXY*H< KKL) /2. 000) 
CEP=GEP*CELX 

KKL, H( KKL) ,EXY,QEP 
CRIT 

HFG , ZFI N ,H (NBOT F ) , TS AT , R PM » QTOT , QT, FANGL 



(6,710) 

(6,720 ) 

(6,730) 

(6,740) 

NR= 1 » NDIV 

(6,7 50 NR, DEL ( NR) , QE(NR) , AMT OT ( NR ) ,TI B ( NR) , TT (N 
1R) ,TE(NR),T6 (NR) 

WRITE (6,760) 

DO 410 NG=1,NDIV,2 

WRITE (6,770) NG, CW ( N G ) , CF ( NG ) , RHOF ( H G 
1 ,R ( NG) ,TBM(NG) ,QINC( NG) 

RETURN 



WRITE 

WRITE 

WRITE 

WRITE 

00 400 

WRITE 



) , UF ( NG ) , EP S( NG ) 




500 

510 
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530 FORMAT (4X,6HCRIT= ,612*5) 

540 FORMAT (2G10.5) 

550 FORMAT (4X,7HFANGL= , E12 . 5 , / ,4X , 5HZ0A= ,E12.5) 

560 FORMAT (15) 

570 FORMAT (4X,5HIFF= ,110) 

580 FORMAT (1615) 

590 FORMAT (2G10.5,4I5) 

600 FORMAT ( ///5X ,4HT I 8 = , I 5 , 10X , 3HTT =, 1 5 , / , 5X ,4HT BM=, I 5 , 10 

1X,3HTE=, I5,/,6X,3HTB=, 15) 

610 FORMAT (// 10X , 17HCRASH , CRASH, CRASH ) 

“ (//5X,4(E12.7,3X) ) 

( 2X, 7HELEMENT, 2 X, 4HQF I N, 1 7X , 7HQ TROUGH ,15X,6HQT0 



620 FORMAT 
630 FORMAT 
1TAL) 

640 FORMAT 
650 FORMAT 



(4X,I5,E12.5,10X,E12.5,10X,E12.5) 

(//,4X, 11HQFIN TOTA L= » E 12 . 5, 10X, 15HQTRCUGH TOTA 
1L= , E12.5 ) 

660 FORMAT (/////, 4X, 5HCL I = , E 12. 5, 5X, 7HCANGL = ,E12.5,/,4X 
1 ,8HRBASE I = ,E12.5,2X,5HR2I= , E12 .5, /, 4X, 8HTH ICKI = ,E12 
2. 5 ,2X » 6HBF I N = , E12. 5 , / ,4X, 5HRPM= , El 2 . 5 , 5X , 5HTSS= ,E12 

3.5, /,4X,6hTINF= , E12 . 5, 4X, 6HH INF = , E 12. 5 , / , 4X , 6HCR I T= 
4.E12.5 ,4X,7HFANGL= , El 2 .5 , / , 4X ,5HZ0A= , E12. 5 , 5X, 5H I FF= 
5 ,110) 

670 FORMAT ( lhl , //2X ,4HB0 A=, G1 2. 5, 5X , 4HZ0 A=, G 12. 5 , 5X, 5HZF I 
1N=,G12.5,5X,5H8FIN=,G12.5,5X ,13HSURFACE AREA=,G12.5 ) 
680 FORMAT ( // 5X , 2HNP , 6X, 1H X, 1 2X , 1HY , 12X , 1HT ) 

690 FORMAT ( /2X , 13 , 3X , 3 ( F 10 .6, 3X ) ) 

700 FORMAT (//2X,2HEL»8X,1HH»11X ,9H EL- LENGTH, 15X , 4HQ-EL ) 
710 FORMAT ( / 2X, 1 2, 3X, E 1 2 . 5, 3X , E 12. 5 , 10X , E 12. 5) 

720 FORMAT ( /2X , 22HCCNV ERGENCE CR ITER I AN= , El 5 .8 ) 

730 FORMAT ( 1H , // , 5X , 4HHFG= , E 12. 5 ,/ , 5X , 1 1 HNO. OF FINS=,E12 

1.5, /,5X,6t-H-0UT = ,E12.5,/,5X,5HTSAT=,E12.5,/, 5X,4HRPM=, 
2E12.5,/ ,5X »6HQ— BGT =,E12.5,/,5X»6HQFIN =, E12 .5 , / , 5X , 11H 
3FALF-ANGLE = , F8. 3) 

740 FORMAT ( 1H0.6X, 1HJ, 4X, 14HFILM TH I CKN ESS , 6X , 8HQ- INCREM , 
16X »8HMASS— TGT »7X,3HTIB»8X,2HTT ,10X ,2HT E , 8X ,2HT B ) 

750 FORMAT ( lh , 4X, 14, 4X , F 12 . 1C, 4X , F 10 .4 , 6X ,F 9. 5 , 6X ,F 5. 1 , 6 
1X*F5.1 ,6X,F5.1,6X,F5.1 ) 

760 FORMAT ( 1HC , 6X , 1H J , 6X, 6HK- WALL ,4X , 6HK-F I LM ,3X , 7HDENS IT 
1 Y, 4X ,9 HV ISC- FILM, 6X,7HEPSI LON, 5X,6HRADIUS, 5X, 3HTBM , 5X , 
25HQ-BCT) 

770 FORMAT ( 1H , 4X, 1 4, 4X , F 7. 3, 4X,F 6. 4, 4X , F 6. 3 ,4X , F9.7 ,4X , F 
19.7, 4X ,F7.5,5X,F5.1, IX, 1P1012.3) 

END 

SUBROUTINE COORD 
IMPLICIT REAL*8 (A-H,0-Z) 

COMMON /GLOB 1/ CLI ,C ANGL , RBA SEI , R2 1 , TH I CKI , BF I N, T Z, TSS 
1 »T INF, HINF, FANGL »ZOA»ZFIN, BOA , QTOT, RPM 
COMMON /PCRO/ X(200) , Y (200 ) , EZERO, 8VIN, THICK, TALFA, APS 
COMMON /APOL/ DOB F , DOTH, KF IN ( 50) , KFF ( 5 0) , I F F , JTC , JLC , J 
1 INT, KT 

DtLH=B VI N/DCBF 
X( 1)=0.0C0 
Y(1)=THICK+8VIN 
N= 1 

CO 20 1=1, IFF 
IC A = KF I N ( I ) 

ICB=KFF< I ) 

CBA=DF LOAT ( I CB-ICA ) 

AN=0. ODO 

OO 10 I 1= ICA, ICB 
X< 1 1 ) = X ( 1 ) +N»AN*DELH*T ALFA/CBA 
Y( 1 1 ) = Y( 1 )-N*DELH 
10 AN=AN+1.0C0 
20 N=N+1 

AN=O.ODO 

ICD= ICB- I CA+1 

DO 50 J=JTC , JLC , JI NT 
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X( J) =X (1 ) 

Y( Jl =( 1. CO C-AN/DOTH)* THICK 
CO 30 JJ=I, ICO 

X( J+JJ)=X(J)+JJ*EZERO/ (2* ( CBA+1 .000 ) ) 

30 Y( J+JJ)=Y( J) 

CO 40 K= 1? KT 

X< J+JJ+K ) = X(J+JJ)+K*APS/ (2.000*KT) 

40 Y( J+J J +K ) = Y ( J ) 

50 AN=AN+ 1 . OCO 
RETURN 
END 

SUBROUTINE FCRMAF 
IMPLICIT REAL*8(A-H,0-Z) 

CIMENS ION E(3),C(3), EA(3,3) „ 

COMMON /GL0B1/ CLI , C ANGL , REAS El , R2 I , T H I CK 1 , 8F IN, T Z , TSS 
I » T IN Ft HINF » F ANGL, ZOA , ZFIN , BOA »QTCT,RPM 
COMMCN / PCRC/ XI200) , Y ( 200 ), EZERO, BVIN, THICK, TALFA, APS 
COMMON /MAFG/ A ( 2C0 , 50 ) , F ( 200 , 1 ) , H (200 ) , TS ( 20 0 ) , TS AT , C 
1K,NEL,NSNP,NBAN, ICOR( 200» 3 ) 

00 20 N= 1 1 NSNP 
F(N. 1)=0.CDC 
DO 10 M A= 1 , N B AN 
10 A(N,MAi=O.CDO 
20 CONTINUE 

CO 70 I EL=1 » NEL 
IA=ICOR( IEL » l) 

IB=ICOR( IEL, 2) 

IC = I CGR ( I EL f 3 ) 

B( 1 ) = Y ( I B ) - Y ( IC) 
fi ( 2 ) =Y (IC)-Y(IA) 

B( 3) =Y ( I A ) — Y ( I B ) 

C(1)=X(IC)-X(IB) 

C(2)=X (IA)-X(IC) 

C ( 3) -X ( I B ) -X ( I A ) 

EL=DSQRT ( C ( 3 )**2+B ( 3 ) **2 ) 

AS=DABS( (B(1)*C(2)-B(2)*C(1))/2»0D0) 

HC=H(I EL) /CK 
DO 60 J= 1 » 3 
JJ=ICGR ( I EL » J ) 

DO 50 K= 1 1 3 
KK = ICOR ( I E L y K ) 

EA(J,K)=(B(J)*B(K)+C ( J )*C ( K) >/(4*AS) 

IF (HC.EC.O.ODO) GO TO 40 
H£L=HC*EL/6.0D0 
IF (J.EQ.2) GO TO 40 

IF (K.EQ.3) GO TO 40 

IF (J.EQ.K) GO TO 30 

EA( J,K)=EA( J,K)+FEL 
GO TO 40 

30 EA(JtK)=EA( J.K)*2*HEL 
40 IF (KK.LT.JJ) GO TO 50 
NW=KK-JJ+1 

A( JJ ,NW) =A( JJ,NW)+EA ( J,K) 

50 CONTINUE 
60 CCNTINUE 

FE=HC*TS( IEL)*EL/2.0D0 
F( IA , 1 ) = F ( IA » 1 ) +FE 
F(IB,1)=F(IB,1)+FE 
70 CONTINUE 
RETURN 
END 

SUBROUTINE BANDEC ( N EQ , MAX 3, NV EC ) 

IMPLICIT REAL*8(A-H,0-Z) 

COMMON /GLOB 1/ CL I , C ANGL , RBA SE I , R2 1 , TH I CKI ,B F I N ,TZ , TSS 
ItTINFtHINFf FANGL, ZOA , ZFIN, BOA, QTCT, RPM 
COMMON /PCRO/ X(200) ,Y(200) , EZERO, BVIN, THICK, TALFA, APS 
COMMON /MAFO/ A ( 200 , 50 ) , F( 200 , 1 ) , H ( 200 ) , TS( 200 ,T SAT ,C 



100 



c 

c 

c 

c 

c 



c 

c 



c 

c 

c 



c 

c 

c 



c 

c 

c 



c 

c 

c 



10 

20 

30 



40 

50 



1K,NEL,NSNP,NBAN,IC0R( 200,3) 

LQGP=N EQ-1 
DO 20 1=1 , LCCP 
MB=I +1 

NB=M INO ( I+MAXB-l.NEQ ) 

DO 20 J=MB , NB 

L= J+2-MB 

D=A( I, L) /A (1,1) 

DO 10 MM=1,NVEC 
F(J,MM)=F(J,MM)-D*F( I , MM ) 
MM=MINO(MAXB-L+i,NEQ-J+l ) 

DO 20 K= 1 , MM 
NN=L+K-1 

A ( J , K) =A( J,K)-D*A( I ,NN) 

DO 30 1= ltNVEC 
F(NEC,I)=F(NEQ,I)/A(NEQ, 1) 

00 5a l=2:,NEQ 
J=NEQ- 1+ 1 

K=MI NO ( NEC- J+l » MAX B ) 

DO 50 MM= 1 ,N VEC 
DO 40 L=2,K 
MB=J+L-1 

F( J, MM J = F ( J » MM )— A ( J » L) *F ( MB , MM) 

F( J , MM ) = F ( J , MM) / A( J » 1) 

RETURN 

END 

SUBROUTINE DPOLRT COMPUTES THE ROOTS OF A REAL 
POLYNOMIAL USING A NEW TON— RA PH SC N ITERATIVE 
TECHNIQUE. 



SUBROUTINE DPOLRT ( X COF, COF, M , ROOTR, ROOT I , IER ) 
IMPLICIT REAL*8 (A-H), REAL*8 (0-Z) 

DIMENS ION XCOF ( l ) , CO F( 1 ) ,RCO TR ( 1 ) .ROOT I ( 1) 



IF IT=0 
N = M 
IER= 0 

IF (XCOF ( N+l ) ) 10,40,10 
10 IF (N) 20,20,60 

SET ERROR CODE TO 1 

20 IER= 1 
30 RETURN 

SET ERROR CODE TO 4 



40 



IER=4 
GO TO 30 



SET ERROR CODE TO 2 



50 IER= 2 
GO TO 30 

60 IF (N-36) 70,70,50 
70 NX=N 

NXX=N+ 1 
N2=l 
KJ 1=N+1 
DO 80 L=1,KJ1 
MT=K J 1— L+ 1 
80 COF( MT )=XCCF (L) 

SET INITIAL VALUES 
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non non ooo non non non 



90 X0=. 00500101 
Y0=0.01000101 

ZERO INITIAL VALUE COUNTER 

IN=0 
100 X=X0 

INCREMENT INITIAL VALUES AND COUNTER 

X0=-10.C*Y0 
Y0=- 10 .0*X 

SET X AND Y TO CURRENT VALUE 

X=X0 
Y=Y0 
IN=IN+1 
GO TO 120 
110 IF IT=1 
XPR=X 
YPR= Y 

EVALUATE POLYNOMIAL AN0 DERIVATIVES 



120 ICT=0 
130 UX=0.0 
UY=0 .0 
V=0.0 
YT=0.0 
XT=1 . 0 
U=CQF( N+ 1 ) 

IF (U) 140,270,140 
140 DO 150 1=1, N 
L=N- 1+ 1 
XT2=X*XT-Y*YT 
YT2=X* YT+ Y*XT 
U=U+CCF (L ) *XT 2 
V=V+CCF(L)*YT2 
F I = I 

UX=UX+FI*XT*COF(L ) 

UY=UY-FI#YT*CCF ( L) 

XT=XT2 
150 YT-YT2 

SUMSQ=UX*LX+UY*UY 
IF ( SUMS Q ) 160, 230, 160 
160 DX = { V*UY-U*UX ) /SUMS Q 
X= X+ D X 

DY=-(U*UY*V*UX I/SUMSQ 
Y=Y+DY 

IF (DABS(0Y)+DABS(DX)-1.0E-05) 210,170,170 

STEP ITERATION COUNTER 

170 ICT=ICT+1 

IF (ICT-5CC) 130,180,180 
180 IF ( IF IT 1 210,190,210 
190 IF (IN-5) 100,200,200 

SET ERROR CODE TO 3 



200 

210 



I ER= 3 
GO TO 30 
DO 220 L= 1 * NXX 
MT=KJ1-L+1 
TEMP=XCGF ( MT ) 
XC0F(MT)=C0F ( L) 
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220 



230 

240 

250 

260 



270 

280 

290 

300 

310 



320 

330 



COF(L)=TEMP 
IT EMP=N 
N=NX 

NX=I TEMP 

IF ( IF IT ) 250, 110, 250 
IF ( IF IT) 240,100,240 
X= XPR 
Y= YPR 
IF I T=0 

IF ( DABS( Y/XJ-1.0E-04) 280,260,260 

ALPHA=X+X 

SUMSQ=X*X+Y*Y 

N=N-2 

GO TO 290 

X= 0. 0 

NX=NX— 1 

NXX=NX X— 1 

Y=O.C 

SUMS Q=0 .0 

ALPHA=X 

N=N-1 

COF (2) =CCF(2)+ALPHA*C0F(1) 

DO 300 L=2,N 

C0F<L+1)=CCF(L+1)+ALPHA*C0F( L J-SUMS3 *COF ( L- 1) 

ROOTI ( N2 ) =Y 

R00TR(N2)=X 

N2=N2+1 

IF (SUMSG) 320,330,320 
Y= -Y 

SUMS G=0. 0 
GO TO 310 
IF (N) 30,30,90 
END 
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